/*
 * Copyright (C) 2014 Nathan Templon
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package mathlib;

import java.io.File;
import java.io.FileWriter;
import java.io.PrintWriter;
import java.text.DecimalFormat;

/**
 *
 * @author Nathan Templon
 */


public class Tester {

	 private static File outputFile;
	 private static FileWriter fw;
	 private static PrintWriter pw;
	 
	// Define given values
	private static final double etaD = 0.97;
	private static final double gammaD = 1.4;
	private static final double etaC = 0.85;
	private static final double gammaC = 1.37;
	private static final double PrB = 1.0;
	private static final double gammaB = 1.35;
	private static final double QR = 4.5E7;
	private static final double etaT = 0.9;
	private static final double gammaT = 1.33;
	private static final double etaN = 0.98;
	private static final double gammaN = 1.36;
	private static final double R = 287; // J/(kg K)
	
	private static final double CpC = (gammaC * R) / (gammaC - 1.0); // J/(kg K)
	private static final double CpB = (gammaB * R) / (gammaB - 1.0); // J/(kg K)
	private static final double CpT = (gammaT * R) / (gammaT - 1.0); // J/(kg K)
	
	private static DecimalFormat format = new DecimalFormat("0.00000E0");
	
    public static void main(String[] args) {
//        openFile();
//        
//        calculateFirstPage("SEA_LEVEL", 101325, 288.16, 0.0, 0.3, 30.0, 1500.0, 1.74844E0);
//        out("\n\n");
//        calculateFirstPage("CRUISE", 24540, 220.02, 0.8, 0.3, 30.0, 951.48113, 1.74844E0);
//        calculateSecondPage(1.95878E-5, 1.93954E-5);
//        
//        closeFile();
    	
    }
    
    public static void calculateFirstPage(String location, double Pa, double Ta, double machInf, double mach02, double PRC, double T04, double A02Cruise) {

        
        out("CpC: " + format.format(CpC));
        out("CpB: " + format.format(CpB));
        out("CpT: " + format.format(CpT));
        
        out("\n" + location + ":\n");
        out("PRESSURE: " + format.format(Pa));
        out("TEMPERATURE: " + format.format(Ta));
        
        out("MACH_INFINITY: " + format.format(machInf));
        
        final double uInf = machInf * Math.sqrt(gammaD * R * Ta); // m / s
        out("U_INFINITY: " + format.format(uInf));
        
        final double T02 = Ta * (1 + 0.5 * (gammaD - 1) * machInf * machInf); // K
        out("T02: " + format.format(T02));
        
        final double P02 = Pa * Math.pow(1 + etaD * ((T02 / Ta) - 1), gammaD / (gammaD - 1)); // Pa
        out("P02: " + format.format(P02));
        
        final double rho02 = P02 / (R * T02); // kg / m^3
        out("Rho02: " + format.format(rho02));
        
        out("MACH_02: " + format.format(mach02));
        
        final double u02 = mach02 * Math.sqrt(gammaD * R * T02); // m / s
        out("U02: " + format.format(u02));
        
        out("PRC: " + format.format(PRC));
        
        final double P03 = PRC * P02; // Pa
        out("P03: " + format.format(P03));
        
        final double T03 = T02 * (1.0 + (1 / etaC) * (Math.pow(PRC, (gammaC - 1) / gammaC) - 1)); // K
        out("T03: " + format.format(T03));
        
        final double P04 = PrB * P03; // Pa
        out("P04: " + format.format(P04));
        
        out("T04: " + format.format(T04));
        
        final double deltaS = CpB * Math.log(T04 / T03) - R * Math.log(P04 / P03); // J / (kg K)
        out("DELTA_S: " + format.format(deltaS));
        
        final double f = (T04 - T03) / ((QR / CpB) - T04);
        out("F: " + format.format(f));
        
        final double T05 = T04 - (1 / (1 + f)) * (CpC / CpT) * (T03 - T02); // K
        out("T05: " + format.format(T05));
        
        final double P05 = P04 * Math.pow(1.0 - (1.0 / etaT) * (1.0 - T05 / T04), gammaT / (gammaT - 1)); // Pa
        out("P05: " + format.format(P05));
        
        final double ue = Math.sqrt(2 * etaN * ((gammaN * R) / (gammaN - 1)) * T05 * (1.0 - Math.pow(Pa / P05, (gammaN - 1) / gammaN))); // m / s
        out("UE: " + format.format(ue));
        
        final double TbyMDotA = (1.0 + f) * ue - uInf; // (N s) / kg
        out("T / mDotA: " + format.format(TbyMDotA));
        
        double mDotA = 0.0; // kg / s
        final double A02; // m^2
        if (location.equals("SEA_LEVEL")) {
            mDotA = 200000.0 / TbyMDotA;
            out("M_DOT_A: " + format.format(mDotA));
            
            A02 = mDotA / (rho02 * u02);
            out("A02: " + format.format(A02));
        }
        else if (location.equals("CRUISE")) {
            A02 = A02Cruise;
            out("A02: " + format.format(A02));
            
            mDotA = rho02 * u02 * A02;
            out("M_DOT_A: " + format.format(mDotA));
            
            final double thrust = TbyMDotA * mDotA; // N
            out("THRUST: " + format.format(thrust));
        }
        
        final double TSFC = f / TbyMDotA;
        out("TSFC: " + format.format(TSFC));
        
        out("");
        
        double T7 = T05 * (1 - etaN * (1 - Math.pow(Pa / P05, (gammaN - 1) / gammaN ) ) ); // K
        out("T7: " + format.format(T7));
        
        double compWork = CpC * (T03 - T02); // J / kg
        out("COMP_WORK: " + format.format(compWork));
        
        double turbWork = CpT * (T05 - T04); // J / kg
        out("TURB_WORK: " + format.format(turbWork));
        
        double keIncrease = 0.5 * (ue * ue - uInf * uInf);
        out("KE_INC: " + format.format(keIncrease));
        
        double etaTherm = 1 - Ta / T03;
        out("ETA_THERM: " + format.format(etaTherm));
        
        double etaTotal = Double.NaN;
        out("ETA_TOTAL: " + format.format(etaTotal));
    }
    
    public static void calculateSecondPage(double TSFC_Takeoff, double TSFC) {
        // Givens
        final double takeOffWeight = 550000.0; // N
        final double finalWeight = 450000.0; // N
        final double CD0 = 0.02;
        final double e = 0.8;
        final double b = 30.0; // m
        final double S = 90.0; // m^2
        final double uClimb = 100.0; // m / s
        final double g = 9.8067; // m / s^2
        final double hCruise = 10500; // m
        final double hMax = 15000; // m
        
        out("\n\nPAGE TWO:\n");
        
        final double AR = (b * b) / S;
        out("ASPECT_RATIO: " + format.format(AR));
        
        final double rhoSL = 1.225; // kg / m^3
        out("RHO_SL: " + format.format(rhoSL));
        
        final double qSL = 0.5 * rhoSL * uClimb * uClimb;
        out("Q_SL: " + format.format(qSL));
        
        final double cLTakeoff = takeOffWeight / (qSL * S);
        out("CL_TAKEOFF: " + format.format(cLTakeoff));
        
        final double cDTakeoff = CD0 + (cLTakeoff * cLTakeoff) / (e * Math.PI * AR);
        out("CD_TAKEOFF: " + format.format(cDTakeoff));
        
        final double dragSL = cDTakeoff * qSL * S; // N
        out("DRAG_SL: " + format.format(dragSL));
        
        final double ROC_SL = (200000 - dragSL) * (uClimb / takeOffWeight); // m/s
        out("ROC_SL: " + format.format(ROC_SL));
        
        final double ROC_Climb = ROC_SL * (1.0 + (hMax - hCruise) / hMax) * 0.5; // m/s
        out("ROC_CLIMB: " + format.format(ROC_Climb));
        
        final double climbWeight = 200000.0 * TSFC_Takeoff * (hCruise / ROC_Climb) * g;
        out("W_CLIMB: " + format.format(climbWeight));
        
        out("");
        
        final double machInf = 0.8;
        out("MACH_INF: " + format.format(machInf));
        
        final double rhoA = .38857; // kg / m^3
        out("RHO_INF: " + format.format(rhoA));
        
        final double Ta = 220.02; // K
        out("TA: " + format.format(Ta));
        
        final double uInf = machInf * Math.sqrt(gammaD * R * Ta); // m / s
        out("U_INF: " + format.format(uInf));
        
        final double initialWeight = takeOffWeight - climbWeight; // N
        out("W_INIT: " + format.format(initialWeight));
        
        final double averageWeight = 0.5 * (initialWeight + finalWeight); // N
        out("W_AVE: " + format.format(averageWeight));
        
        final double CLclimb = averageWeight / (0.5 * rhoA * uInf * uInf * S);
        out("CL_CLIMB: " + format.format(CLclimb));
        
        final double CDclimb = CD0 + (CLclimb * CLclimb) / (e * Math.PI * AR);
        out("CD_CLIMB: " + format.format(CDclimb));
        
        final double LDclimb = CLclimb / CDclimb;
        out("L/D CLIMB: " + format.format(LDclimb));
        
        final double TRclimb = averageWeight / LDclimb; // N
        out("TR_CLIMB: " + format.format(TRclimb));
        
        out("");
        
        final double Smax = (uInf / (TSFC * g)) * LDclimb * Math.log(initialWeight / finalWeight); // m
        out("S_MAX: " + format.format(Smax));
        
        final double maintFact = 1.0;
        out("MAINT_FACTOR: " + format.format(maintFact));
        
        final double costPerHr = 320.0 + 500.0 * maintFact; // $ / hr
        out("COST_PER_HOUR: " + format.format(costPerHr));
        
        final double Sroute = 2.7E6;
        out("S_ROUTE: " + format.format(Sroute));
        
        final double flightTime = (((hCruise / ROC_Climb) * 2) + (Sroute / uInf)) / 3600.0; // hr
        out("FLIGHT_TIME: " + format.format(flightTime));
        
        final double cruiseWeight = initialWeight * (1 - (1 / Math.exp((Sroute * TSFC * g) / (uInf * LDclimb)))); // N
        out("W_CRUISE: " + format.format(cruiseWeight));
        
        final double fuelCost = 3 * (climbWeight + cruiseWeight) * (1 / 30.426); // $
        out("FUEL_COST: " + format.format(fuelCost));
        
        final double seatKmCost = (costPerHr * flightTime + fuelCost + 3000) / ((Sroute / 1E3) * 100);
        out("SEAT_KM_COST: " + format.format(seatKmCost));
    }
    
    public static void openFile() {
        try {
        outputFile = new File("D:\\Dropbox\\output.txt");
        fw = new FileWriter(outputFile);
        pw = new PrintWriter(fw);
        }
        catch(Exception e) {
            System.out.println(e);
        }
    }
    
    public static void closeFile() {
        try {
        pw.close();
        fw.close();
        }
        catch(Exception e) {
            System.out.println(e);
        }
    }
    
    public static void out(Object o) {
        pw.println(o.toString());
    }
}
